import numpy as np
import matplotlib.pyplot as plt

a1 = np.loadtxt('Cl401.txt')
a2 = np.loadtxt('S193.txt')
a3 = np.loadtxt('S313.txt')
a4 = np.loadtxt('S329.txt')
a5 = np.loadtxt('S349.txt')

# S329-S349-Cl401
# S329-S349-S313
# S329-S349-S193
# 7 bonds in total

#print(len(a1))

def triangle_area(a,b,c):
    s = (a + b + c) / 2
    area = (s*(s-a)*(s-b)*(s-c))** 0.5
    return area


# Area of a triangle = (s*(s-a)*(s-b)*(s-c))** 0.5
# s = (a + b + c) / 2

area1 = []
area2 = []
area3 = []

for i in range(len(a1)):
    l1 = np.linalg.norm(a1[i]-a4[i])  # Cl401-S329
    l2 = np.linalg.norm(a1[i]-a5[i])  # Cl401-S349 
    l3 = np.linalg.norm(a4[i]-a5[i])  # S329-S349
    l4 = np.linalg.norm(a3[i]-a4[i])  # S313-S329
    l5 = np.linalg.norm(a3[i]-a5[i])  # S313-S349 
    l6 = np.linalg.norm(a2[i]-a4[i])  # S193-S329
    l7 = np.linalg.norm(a2[i]-a5[i])  # S193-S349

    area1.append(triangle_area(l1,l2,l3))
    area2.append(triangle_area(l3,l4,l5))
    area3.append(triangle_area(l3,l6,l7))

area1 = np.array(area1)
area2 = np.array(area2)
area3 = np.array(area3)

np.savetxt('area1.txt',area1,fmt='%5.5f')
np.savetxt('area2.txt',area2,fmt='%5.5f')
np.savetxt('area3.txt',area3,fmt='%5.5f')

t = np.linspace(0,len(a1)/100,len(a1))
np.savetxt('t.txt',t,fmt='%5.5f')

plt.plot(t,area1)
plt.plot(t,area2)
plt.plot(t,area3)
plt.xlabel('t (ps)')
plt.ylabel(r'Area ($\AA^2$)')
plt.ylim(4,10)
plt.legend(['S329-S349-Cl401 T5','S329-S349-S313 T2','S329-S349-S193 T5'],frameon=False)
plt.savefig('area.png',bbox_inches='tight')
plt.show()
